Ab initio theory of gate induced gaps in graphene bilayers 
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We study the gate-voltage induced gap that occurs in graphene bilayers using ab initio density 
functional theory. Our calculations confirm the qualitative picture suggested by phenomenological 
tight-binding and continuum models. We discuss enhanced screening of the external interlayer 
potential at small gate voltages, which is more pronounced in the ab initio calculations, and quantify 
the role of crystalline inhomogeneity using a tight-binding model self-consistent Hartree calculation. 



I. INTRODUCTION 

Recently, ultrathin graphite films including 
monolayers^ 3 and bilayers-'- have attracted con- 
siderable attention because of their novel properties. In 
single-layer graphene, the A sublattice to B sublattice 
hopping amplitude vanishes at two incquivalent points 
K and K' on the edge of the honeycomb lattice Brillouin 
zone (BZ); away from these points, the hopping ampli- 
tude grows linearly with the wave vector and has a phase 
which winds along with the orientation of the wave 
vector measured from the high-symmetry points. The 
band structure of an isolated graphene layer is therefore 
described at low energies by a two-dimensional massless 
Dirac equation with linear dispersion; this property 
gives rise to a half integer quantum Hall effecl£i£ and 
to a quantized spin Hall effect^ and dominates the 
low-energy physics. 
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gate bias have been studied experimentally using angle- 
resolved photoemission spectroscopy 5 and Shubnikov-de 
Haas analysis^ of magnetotransport. This unique prop- 
erty of bilayer graphene has created considerable interest 
in part because it suggests the possibility of switching 
the conductance of a graphene bilayer channel over a 
wide range at a speed which is limited by gate-voltage 
switching, as illustrated schematically in Fig. [3] 
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FIG. 2: (Color online) Schematic illustration of a circuit with 
a bilayer graphene channel sensitive to an external gate volt- 
age. The graphene channel is separated from the front and 
back gates by a Si02 layer. Fhe channel resistance change will 
be rapid and large when the graphene channel is undoped and 
isolated from the gate electrodes, as illustrated here. In this 
case, the total charge density in the bilayer system is fixed 
and the chemical potential lies in the gap opened by the gate 
voltage. This geometry could also be used to capacitively 
probe the correlation physics of the isolated bilayer system, 
as discussed in the text. 



FIG. 1: (Color online) Structure of a graphene bilayer with 
honeycomb lattice constant a = 2.46 A and interlayer sepa- 
ration d = 3.35 A. 

In bilayer graphene, the Bernal (A-B) stacking illus- 
trated in Fig. [T] modifies this electronic structure in an 
interesting way.— ^ At K and K' , the states localized at 
the A and B sites, are repelled from zero energy by in- 
terlayer tunneling; only states localized at A and B are 
present at zero energy. When tunneling is included, the A 
to B hopping is a second-order process via a virtual bond- 
ing or antibonding state at A and B. The chirality of 
the low-energy bands is therefore doubled. Most intrigu- 
ingly, an external potential which induces a difference be- 
tween the A and B site energies will open up a gap 11 ! 12 
in the spectrum. Band gaps controlled by applying a 



In this paper, we report on an ab initio density func- 
tional theory (DFT) study of the influence of an external 
potential difference between the layers on the electronic 
structure of a graphene bilayer. We compare our results 
with the phenomenological tight-binding and continuum 
model Schrodinger-Poisson calculations used in previ- 
ous theoretical analyses.— DFT predicts, in agreement 
with these works, that the external potential difference 
is strongly screened with a maximum energy gap value 
of ^0.3 eV. There are, however, quantitative differences. 
In particular, the enhanced screening which occurs for 
weak external potentials is stronger in the DFT calcula- 
tions. In an effort to improve the quantitative agreement, 
we have estimated the influence of crystalline inhomo- 
geneity in a tight-binding model self-consistent Hartree 
calculation. This effect strengthens intralayer Coulomb 
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interactions because the charge is spatially bunched, and 
therefore increases screening in a Hartree calculation, but 
does not fully account for differences between the two cal- 
culations. 



II. AB INITIO DENSITY FUNCTIONAL 
THEORY CALCULATIONS 

We have performed ah initio DFT calculations^ for 
an isolated graphene bilayer under a perpendicular ex- 
ternal electric field using an all-electron linearized aug- 
mented plane wave plus local-orbital method incorpo- 
rated in WIEN2K. 15 We used the generalized-gradient 
approximation^ for the exchange and correlation poten- 
tial. 



A. External electric fields 

To investigate the influence of an external electric field 
on a graphene bilayer, a periodic zigzag potential was ap- 
plied along the z direction, perpendicular to the graphene 
planes, in a supercell. 17 The bilayer was placed at the cen- 
ter of the constant external electric field region and the 
size of the supercell was set to a large value (~16 A) to 
minimize the interaction between bilayers in neighboring 
supercells. In order to resolve the small gaps produced by 
small external fields we performed BZ sums using a rel- 
atively large number of k points (~800) per irreducible 
wedge (5000 k points in the whole BZ). Total energies 
were convergent to within 0.0001 Ry. 




FIG. 3: (Color online) An averaged Coulomb potential of a 
cross section vs z for total external potential U ex t = eV and 
Uext = 1 eV. Here, zo is the superlattice period and the cross 
section was chosen to include equal number of atoms in each 
layer. 



Figure [3] shows the Coulomb potential relative to the E 



Fermi energy, laterally averaged along a line in the x-y 
plane that includes an equal number of atoms in each 
layer, as a function of the z coordinate. The potential 
includes the Hartree electron-electron potential and the 
electron-ion interaction but not the external electric field 
potential or the exchange-correlation potential. The bi- 
layer is centered around z/zq = 0.25, where zq is the 
superlattice period. (In the discussion below, we define 
the external potential energy U ex t as U ex t = eE zext d, 
where E z . ext is an external electric field along the z direc- 
tion and d is the interlayer separation of bilayer graphene 
which we take to be 3.35 A.) In the absence of an elec- 
tric field (dashed line), the Coulomb potential is flat in 
the vacuum region and the energy difference between the 
vacuum and the Fermi energy gives estimates of the work 
function of bilayer graphene to be ^4.3 eV. In the pres- 
ence of an electric field (solid line), charge transfer be- 
tween the layers induces a potential which cancels the 
external potential in the vacuum region. The difference 
between the Coulomb energies of the two layers in the 
presence of an external electric field is closely related to 
the gate- voltage induced energy gap. 



B. Energy bands 

Figure 0Ja) shows the DFT energy band structure 
of bilayer graphene in the absence of an applied exter- 
nal electric field. When U ex t = eV, the low-energy 
band dispersion is nearly parabolic at two inequivalent 
corners, K and K' , of the hexagonal BZ, as predicted 
by the 7r-orbital tight-binding and continuum model 
phenomenologies . 1 1 ' 12 The valence and conduction bands 
meet at the Fermi level. 

In the absence of an external electric field, bilayer 
graphene, like single-layer graphene, is a zero-gap semi- 
conductor. At finite U ex t, however, the low-energy bands 
near the K or K' point split, as explained in the In- 
troduction. Therefore, gated graphene bilayer systems 
are gate- voltage tunable narrow gap semiconductors [Fig. 
HJb)]. This property is unique, to our knowledge. It is 
worth noting that in the presence of an external electric 
field, the true energy gap does not occur at the K or K' 
point but slightly away from it. The low-energy spec- 
trum develops a Mexican hat structure as the strength 
of the external electric field increases. This property is 
also captured by phenomenological models of graphene 
bilayers i^ 2 - 



Evolution of tight-binding model parameters 

With Uext 



Figure [5] illustrates DFT predictions for the evolution 
of tight-binding parameters with the applied external po- 
tential. The tight-binding model expression for the four 
low-energy band eig envalues at the K and K 1 points is 



K/K' 



±U/2, ±\/ji + U 2 /4P- where U is the inter- 
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FIG. 4: (Color online) (a) Bilayer graphene band structure in 
the absence of an external electric field, (b) Bilayer graphene 
band structure near the K point for U ex t=0, 0.5, and 1 eV. 



layer energy difference and 71 is the interlayer tunneling 
amplitude. (As we discuss below, this expression should, 
strictly speaking, be slightly modified in the presence of 
an external potential, but it still provides a convenient 
way of characterizing DFT predictions for the low-energy 
bands.) The values of U and 71 plotted in Fig. \5\ repre- 
sent this interpretation of the four lowest-energy DFT 
eigenvalues and clearly reflect substantial screening of 
the external interlayer potential by the Hartree potential 
plotted in Fig. [3J The interlayer coupling 71 increases 
monotonically as the external potential increases. The 
rate of increase is, however, ten times smaller than es- 
timated in a recent experimental study^ of a doped bi- 
layer systems, possibly suggesting significant differences 
between doped and undoped systems. The intralayer 
nearest-neighbor 7r-electron hopping amplitude 70 and 
the interlayer A-B coupling 73 were fitted to reproduce 
the band dispersion around the K/K' points at low en- 
ergies. We find 70 ~ 2.6 eV and 73 » 0.3 eV, nearly in- 
dependent of the external electric field. This value for 70 
corresponds to an in-plane velocity v = ^Jp- « 8.4 x 10 5 






f 




t 

i 
/ 

/ ' 




i 

f 


f 



/ 









0.342 



0.341 

i 
O.340 



0.0 0.5 1.0 

0.4 



1.5 0.0 0.5 1.0 1.5 



0.339 



0.3 




©•© ab tnttio 

— - tignt-bindirtg. y a =0»V 

— tight-binding, v 3 =0.3 eV 



0.4 D.6 

U (sV) 



0.8 



1.0 



FIG. 5: (Color online) (a) Evolution of the graphene bilayer 
screened on-site energy difference U, extracted from the ab 
initio DFT bands as explained in the text, with the external 
potential U ex t- The external potential is strongly screened, 
(b) Evolution of the interlayer tunneling amplitude 71 with 
Uext- (c) Comparison of band gap as a function of the on-site 
energy difference U obtained from the ab initio DFT calcu- 
lations (open circles) with the tight-binding result for 73 = 
(dashed line) and 73 = 0.3 eV (solid line). 



m/s, where the lattice constant a — 2.46 A. 

Figure [5jc) compares the relationship between the on- 
site energy difference U extracted from the DFT calcu- 
lations and the energy gap with the corresponding rela- 
tionship in the tight-binding model. Note that the gap 
does not increase indefinitely with U but saturates at 
^0.3 eV due to the Mexican hat structure shown in the 
bands illustrated in FigJU For 73 = 0, we can estimate 
the approximate energy gap from the low-energy ap- 
proximation of the tight-binding model given by E gap w 
|f/|7i/V7i + U 2 , where E gap approaches 71 ~ 0.34 eV 



as U increases^ For 73 
duced from that of 73 = 



s 0.3 eV, however, E gap is re- 
and matches well with the 
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DFT results. A nonzero value for 73 has a noticeable 
quantitative influence on the bands. This agreement con- 
firms (unsurprisingly) that the tight-binding model cap- 
tures the character of the low-energy bands in bilayer 
graphene. The most interesting physics is in the relation- 
ship between U and U e xt, which we now examine more 
closely. 



III. SCREENING THEORIES 
A. Continuum Hartree potential models 

The screening of the external potential has been ex- 
amined previously for both doped and undoped bilay- 
ers using phenomenological approaches combined with 
the Poisson equation*^ This type of analysis provides 
a good reference point for interpreting the DFT results 
so we start with a discussion of this picture. Consider 
a graphene bilayer with an interlayer separation d un- 
der an external electric field E z ^ ext along the z direction. 
Neglecting the finite thickness and crystalline inhomo- 
geneity of the graphene layers, and screening external to 
the bilayer, the Poisson equation is 
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V • E = 4?r(-e) K<5(z) + n 2 5(z - d)} , 



(1) 



where n\ and n 2 are the net charge densities on the bot- 
tom and top layers, respectively. If the bilayer is placed 
on a gate dielectric such as silicon dioxide (Si02) and a 
voltage is applied between a gate and the bilayer, an ex- 
cess charge carrier density n = m + n 2 is supplied to the 
bilayer graphene and redistributed between the top and 
bottom layers due to an external electric field. 

In order to compare with our DFT calculations, we 
focus here on the isolated bilayer case illustrated in FigJ^l 
in which the total excess density n — n\ + n 2 = 0. Let 
us define Sn — n 2 = —rti. From Eq. (fTJ, we obtain the 
screened electric field E z between the graphene sheets of 
the bilayer to be 



E, — E, 



AiteSn. 



(2) 



Adding the corresponding Hartree potential to the exter- 
nal potential, we obtain the screened interlayer potential 
difference as 



U = U,, 



4:Tre 2 dSn, 



(3) 



where U = e,E z d and U ex t = eE z , ext d. 

To estimate the relationship between U and U ex t, we 
need only a theory for the dependence of Sn on U. In the 
7r-orbital tight-binding model, Sn is given by the follow- 
ing integral over the BZ: 

d 2 k 



Sn= K 2 L 



r (^(fc)|-%(fc)), (4) 



where \tpi(k)) is a band eigenstate in the presence of U, 
a z = diag(l, 1, — 1, — 1) in the (top,bottom)x (A, B) ba- 
sis, and the index i runs over all occupied states. The 
factor of 2 was included to account for spin degeneracy. 
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FIG. 6: (Color online) The ratio of the external electric po- 
tential Ue X t to the interlayer energy difference inferred from 
the ab initio DFT calculation compared with the value of 
the same ratio in tight-binding model self-consistent Hartree 
calculations, both with and without crystalline inhomogene- 
ity corrections. The tight-binding model calculations used 
7o=2.6 eV for the intralayer tunneling amplitude, 71=0. 34 
eV for the interlayer tunneling amplitude, and 73=0.3 eV for 
the interlayer A-B coupling. 



Figure [5] compares the screening ratio U ex t/U obtained 
from the ab initio DFT calculations with the screening 
ratio from 7r-orbital tight-binding model self-consistent 
Hartree calculations with and without corrections that 
account for the crystalline inhomogeneity within each 
layer as explained later. The agreement between the 
three different approaches is generally good especially at 
large potentials. 

Note that as U approaches zero, U ext /U increases in all 
approximations. This property reflects increased screen- 
ing as the gap decreases and is explained most succinctly 
using the two-band continuum models for the lowest- 
energy bands: 



H, 



eff 
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~2 



1 

-1 



2m 



(7ft) 2 
7T 2 



(5) 



p x -p y 



where it = p x + ip y , m = ^2, 

and er are 2x2 Pauli matrices describing the top and bot- 
tom layer low-energy sites. This Hamiltonian has simple 
spectra e± — ±|a| with cigcnfunctions given by 



l p -i<t>/i 



sin ■ 



itf>/2 
72 



(6) 
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where tan 9 



Sn = 4 



and tan < 



d 2 p 



It follows that 



(~,p\^-\-,p) 



|p|<Pc (2^r 2 



= T2 / pdp cos 6 (p) 
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x c + yxl + l 



(7) 



where x c — We have inserted a factor of 4 in this 

continuum model calculation to account for both spin(| 
and |) and valley(i4T and K') degeneracies. The integral 
over wave vector was cut off at the radius p c ~ ^/2rnr/i 
beyond which the continuum model fails. 
Inserting Eq. (JT]) in Eq. ([3]), we obtain 
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where as = h 2 /m e e 2 is the Bohr radius and m e is the 
bare electron mass. For small U, x c is large and this 
simplifies to 



U,, 



u 



—) f^L lln 



mf7 



(9) 



A related observation concerning the logarithmic diver- 
gence of the screening ratio at small gate voltages was 
made previously by McCann^ All three of our calcu- 
lations exhibit this increased screening at weak external 
potentials, with the largest upturn in the ab initio calcu- 
lations. 



B. Lattice Hartree potential models 

We now turn our attention to one important contribu- 
tion to discrepancies between the ab initio DFT results 
and the predictions of self-consistent Hartree models sim- 
ilar to those described above, the role of crystalline in- 
homogeneity in bilayer and single-layer graphcnc electro- 
statics. We consider a general two-body interaction term 
V, 



v = l 



, } J2 (WIVlAxAa) 4, c A , c A2 c Al , (10) 

X 1 ,A 2 , Ai , A2 



where c x and c\ are creation and annihilation opera- 
tors for a state A. To capture the main consequences of 
crystalline inhomogeneity, we assume that the 7r-orbital 
Bloch states with crystal momentum k can be written as 
a linear combination of atomic orbitals, 
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v ]T e ikR cl> x (x-R-T X ) 



(ii) 



where is an atomiclike it orbital, R is a lattice vector, 
t\ is the displacement of the sites in a unit cell with 
respect to the lattice vector, and TV is the number of 
lattice sites. If we assume that the overlap of 0A-orbitals 
centered on different sites can be neglected and ignore the 
z direction spread of the graphene sheets, the interaction 
Hamiltonian simplifies to 



V 



^ E E V Xu x 2 (q) cl i+qM , 
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fei,fe 2 ,q Ai,A 2 

where Q is the area of the two-dimensional plane 



fc2-q,A 2 Cfe 2,A2 C fcl,Ai ! 

(12) 



i,A 2 (g) = J dx 



and 



wx(q) 



dx e 



-iqx 



-iq-x 



Vx 1 ,x 2 (x), 
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(13) 
(14) 



(15) 



Note that the labels fci and &2 are restricted to the BZ, 
while q runs over the two-dimensional plane. In Eq. 
(fl"4|) . V\ lt \ 2 {x) = e 2 /\x\ when Ai and A2 refer to sites 
in the same layer and V\ lt \ 2 {x) — e 2 / \/\x\ 2 + d 2 when 
Ai and A2 refer to layers separated by d. It follows that 
V\ 1 .\ 2 (q) — 2ne 2 /\q\ for labels in the same layer and 
Vx lt x 2 (q) — 27re 2 exp(— |q|ci)/|q| for labels in different 
layers. Since the total charge of the bilayer is fixed in 
our calculations, only the differences between the various 
Vx lt x 2 (q) values are relevant. For explicit calculations, 
we have used a Gaussian form factor w\(q) = e~' 9 ' r °^ 2 
corresponding to |0a(^)|" where ro ~ 0.48 

A was obtained by fitting to the DFT valence orbitals. 

This two-body Hamiltonian can be used to account for 
crystalline inhomogeneity in a graphcnc bilayer system 
with arbitrary electronic correlations. To compare with 
the ab initio DFT calculations, we consider interactions 
in a mean-field Hartree approximation in which the in- 
teraction contribution to the single-particle Hamiltonian 



E 



e (i V a. 



(16) 



where a and a denote layer and sublattice degrees of 
freedom. Here, 



e (H) - v y , ,n > 

c acr / v v a<7,a' a' ' b a> 



(17) 



where n aa = ^ ^ fe (c^^cu^ including spin degener- 
acy and 



(18) 
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with G a triangular lattice reciprocal-lattice vector. 

As explained in the Introduction, interlayer tunneling 
in graphene leads to high-energy bands which favor the 
A-B sites and low-energy bands that favor the A-B sites. 
Since the low-energy bands respond most strongly to the 
external potential, we can expect that the charge trans- 
fer occurs more strongly on the A-B sites, and that the 
screening potential should be larger on these sites. In- 
stead of a single-interlayer Hartree screening potential, 
two Hartree potentials for low and high bands must be 
calculated separately: 



(H) 
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SH) _ JH) 
A > 

(H) 
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(19) 



When only the G = 
lattice vector sum, 



term is retained in the reciprocal- 
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Am 



2ne 2 d (n^ + rig — tia — n B ), (20) 



and Eq. ([3]) is recovered. It turns out that the sum 
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FIG. 7: Splitting of Hartree potentials (e, (H) - e[ H, )/(e l 
eji ) as a function of (a) the external electric potential U ex t 
and (b) the corresponding density inhomogeneity (Am — 
Ani,) / (An; + An^) in the lattice Hartree potential model. 

over reciprocal-lattice vectors can be truncated with good 
accuracy at the first shell. Noting that e~' G ' d <C 1, we 
find for the crystalline inhomogeneity corrections 



c\ H1) sa 2Tre 2 da(G)(6Ani 
« 2ne 2 da(G)(6An h 

is - ha, An h = n A - 



■ 3An h ), 
-3Ani), 
riB, and ct(G) 



(21) 



where An; = n 

e -\ G \ 2r o /\G\d fa 0.0136. Thus, the inhomogeneity effect 
results in more screening as expected, but as indicated by 
the black dashed line in Fig. [HI it is not able to account 
for the largest part of the discrepancy between DFT and 
model results. 



As the external electric potential is decreased, the dif- 
ference between low-energy and high-energy site occu- 
pancies is increased. The difference in Hartree potentials 
rises correspondingly, as illustrated in Fig. EJa). From 
Eq. (|2ip . we can estimate the relation between the split- 
ting of the Hartree potentials and the density inhomo- 
geneity: 



(H) _ AH) 



9 ,,-An, 
2 a(G W 



An h 



Ann 



(22) 



where the coefficient §a(G) is given by ~0.0612 [Fig. 
Mb)}- 



IV. DISCUSSION 

Our DFT calculations of external potential induced 
gaps in the electronic structure of graphene bilayers con- 
firm the simple picture provided by phenomenological 
tight-binding models. The ab initio calculations include 
a number of effects not contained in the model calcu- 
lations. For example, the occupied a orbitals within 
each graphene plane, which are neglected in the 7r-orbital 
tight-binding model, will be slightly polarized by the ex- 
ternal electric field and contribute to screening. In the 
DFT calculations, not only Hartree potentials but also 
exchange-correlation potentials will be altered by an ex- 
ternal electric field and influence the screening process. 
Since the exchange potential is attractive, its contribu- 
tion to the total potential will lower energies in a layer 
more as the density is increased. The exchange potential 
therefore makes a negative contribution to the screening 
ratio. The quantitative discrepancies between the DFT 
and phenomenological model reflect the combination of 
these and other additional effects contained in the DFT 
calculations, and strong sensitivity to intralayer and in- 
terlayer tunneling amplitudes which may not be evalu- 
ated with perfect accuracy by DFT. We also note that the 
low-energy eigenstates in bilayer graphene are coherent 
combinations of amplitudes on both layers, which implies 
that interlayer exchange interactions will be substantial. 
This kind of effect is absent in the exchange-correlation 
potentials commonly used in DFT. Indeed, it is entirely 
possible that DFT calculations do not predict accurate 
values for the screening ratio. We believe that there is 
strong motivation for capacitive studies of the interlayer 
screening properties of graphene bilayers using an exper- 
imental arrangement similar to that in Fig. fj] 

In summary, we have used ab initio density func- 
tional theory calculations to study the gate-voltage tun- 
able gap in the electronic structure of bilayer graphene. 
The electric-field dependence of the on-site energy dif- 
ference and the interlayer tunneling amplitude were ex- 
tracted from the DFT calculation results by fitting to 
tight-binding model expressions for high-symmetry point 
graphene bilayer band eigenvalues. The screening effect 
seen in the DFT calculations can be explained by a tight- 
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binding model self-consistent Hartree method including 
crystalline inhomogeneity corrections, although the DFT 
screening is stronger especially for weak external poten- 
tials. 
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